	********************************************************************************
	**** Date: 5.03.2019														****
	**** Update: 2.09.2022 														****
	**** Author: JW josephgwright@gmail.com 									****	
	** NOTE: This program has been executed in Stata 17.0						****
	***************** **************************************************************

	***********************************************************************************
	**** Additional modules and packages required for executing the following code ****
	***********************************************************************************
		** btscs.ado
		** ciplot.ado
		** estout.ado
		** ivreg2.ado	
		** krls.ado
		** listtex.ado
		** reghdfe.ado
		** sutex.ado
		** xtivreg2.ado
		** crossfold.ado
		** cvauroc.ado

	******************************
	**** Set directory, seed *****
	******************************
		global dir = "C:\Users\jgw12\Dropbox\Research\Pers-NAVCO\Pers-RENAVCO\CSW-BJPS-reproduction"
		cd "$dir"
		
		capture log close
		log using Pers-Protest-Analysis.log, replace
		set scheme plotplain
		set more off 
		set matsize 1000
		global seed ="984353"

			
 ************************************************************************************
 ************************************   Analysis   **********************************
 ************************************************************************************
		use temp-fe,clear
		global cvar = "coldwar lxyrs lt lnregion"
		qui centile xpers2 if sample==1,centile(50)
		di r(c_1)
		gen treat  = xpers2>=-.02472032 if xpers2~=.
		tab treat if sample==1
		egen m_treat = mean(treat) if sample==1,by(gwf_caseid)
		egen min=min(year) if treat==1,by(gwf_caseid)
		gen first_treatment =year==min if year~=.
		sort cowcode year
		keep if year~=. & xpers2~=. & gdem~=.
		save temp-id,replace
		
		* Show within equivalency in OLS *
			xtset gwf_caseid year
				* unit means *
			qui reg xonset lxyrs lt lnregion treat coldwar m_lxyrs m_lt m_lnregion m_coldwar m_treat m_xonset, vce(cluster gwf_caseid)
			lincom treat
				* unit intercepts *
			qui reg xonset lxyrs lt lnregion treat coldwar i.gwf_caseid, vce(cluster gwf_caseid)
			lincom treat		
				* fixed effects *
			qui xtreg xonset lxyrs lt lnregion treat coldwar,fe vce(cluster gwf_caseid)
			lincom treat
				* absorbed fixed effects *
			qui reghdfe xonset lxyrs lt lnregion treat coldwar,a(gwf_caseid) vce(cluster gwf_caseid)
			lincom treat
				* unit means without \bar{Y} *
			qui reg xonset lxyrs lt lnregion treat coldwar m_lxyrs m_lt m_lnregion m_coldwar m_treat, vce(cluster gwf_caseid)
			lincom treat
			* Nonlinear link *
			qui probit xonset lxyrs lt lnregion treat coldwar m_lxyrs m_lt m_lnregion m_coldwar m_treat, vce(cluster gwf_caseid)
			margins,dydx(treat)
			lincom treat
			qui xtprobit xonset lxyrs lt lnregion treat coldwar m_lxyrs m_lt m_lnregion m_coldwar m_treat, vce(cluster gwf_caseid)
			margins, dydx(treat) predict(pu0)
			lincom treat
 			
			 * ECM LRM *
			xtset gwf_caseid year
			ivreg  xonset (d.xonset=l.xonset) d.lxyrs lxyrs d.lt lt d.lnregion lnregion d.treat treat, cluster(gwf_leaderid)
			ivreg  xonset (d.xonset=l.xonset) d.lxyrs lxyrs d.lt lt d.lnregion lnregion d.xpers2 xpers2, cluster(gwf_leaderid)
			
			* LDV-FE * 
			xtset cowcode year
			gen l1mean5 = l.mean5
			alpha l1v2cademmob l1mean5 lag_xongoing  l1_nav13_ongoing l1_nav21_ongoing, std item gen(l1protest)
			gen l2protest = l.l1protest
			reghdfe xonset lxyrs l1protest l2protest lt lnregion treat,a(gwf_caseid year)vce(cluster gwf_leaderid)
			reghdfe xonset lxyrs l1protest l2protest lt lnregion xpers2,a(gwf_caseid year)vce(cluster gwf_leaderid)
			
			* FD *
			xtset gwf_leaderid year
			xi:xtivreg2 xonset lxyrs lt lnregion treat,fd cluster(gwf_leaderid)  
			xi:xtivreg2 xonset lxyrs lt lnregion xpers2,fd cluster(gwf_leaderid)  
			
			* Leader-FE  *
			reghdfe xonset lxyrs lt lnregion treat,a(gwf_leaderid year)cluster(gwf_leaderid)
			reghdfe xonset lxyrs lt lnregion xpers2,a(gwf_leaderid year)cluster(gwf_leaderid)
			

		* Base model *
		use temp-fe, clear
		xtset gwf_caseid year
		probit xonset $cvar xpers2,vce(cluster gwf_caseid)
		margins,dydx(xpers2 lt lnregion)nose
		est store b0
		xtprobit xonset $cvar xpers2,vce(cluster gwf_caseid)
		margins,dydx(xpers2 lt lnregion)nose
		est store b1
 		xtprobit xonset $cvar xpers2 m_coldwar m_lxyrs m_lt m_xpers2 m_lnregion,vce(cluster gwf_caseid)
		margins,dydx(xpers2 lt lnregion)nose
 		est store b2
		xtprobit  xonset lxyrs lt lnregion xpers2 m_lxyrs m_lt m_xpers2 m_lnregion y_lxyrs y_lt y_xpers2 y_lnregion,vce(cluster gwf_caseid)
		margins,dydx(xpers2 lt lnregion)nose
		est store b3
		xtprobit xonset lxyrs lt lnregion xpers2 m_lxyrs m_lt m_xpers2 m_lnregion y_lxyrs y_lt y_xpers2 y_lnregion mXy_lxyrs mXy_lt mXy_xpers2 mXy_lnregion if sample==1,vce(cluster gwf_caseid)
		margins,dydx(xpers2 lt lnregion)nose 
		est store b4
		
		* Check with IRT-2PL pers instead of generalized *
		qui sum irtpers if sample==1
		replace irtpers=(irtpers-r(mean))/r(sd)
		sum irtpers
		qui egen m_irtpers =mean(irtpers2) if sample==1,by(gwf_caseid)
		qui egen y_irtpers =mean(irtpers2) if sample==1,by(year)
		xtprobit xonset $cvar irtpers2 m_coldwar m_lxyrs m_lt m_irtpers m_lnregion,vce(cluster gwf_caseid)
		margins,dydx(irtpers)nose
		xtprobit xonset lxyrs lt lnregion irtpers2 m_lxyrs m_lt m_irtpers m_lnregion y_lxyrs y_lt y_irtpers y_lnregion,vce(cluster gwf_caseid)
		margins,dydx(irtpers)nose

		
		* Various measures of personalism *
		xtset gwf_caseid year
		probit xonset $cvar gwf_pers,vce(cluster gwf_caseid)
		est store p1
		probit xonset $cvar xpers,vce(cluster gwf_caseid)
		margins,dydx(xpers lt)
		est store p2
		xtprobit xonset  $cvar xpers,vce(cluster gwf_caseid)
		margins,dydx(xpers lt)
		est store p3
		xtprobit xonset $cvar xpers1,vce(cluster gwf_caseid)
		margins,dydx(xpers1 lt)
		est store p4	
		xtprobit xonset $cvar xpers2,vce(cluster gwf_caseid)
		margins,dydx(xpers2 lt lnregion)
		est store p5			
		xtprobit xonset $cvar xpers1 xpers2,vce(cluster gwf_caseid)
		est store p6
		xtprobit xonset $cvar gwf_pers,vce(cluster gwf_caseid)
		
		* Random intercept *
		xtset gwf_caseid year
		melogit xonset $cvar xpers2||gwf_caseid:,vce(cluster gwf_caseid)	
		* Random slope *
		melogit xonset $cvar xpers2||gwf_caseid:xpers2,vce(cluster gwf_caseid)
		melogit xonset $cvar xpers2||gwf_caseid:lt lxyrs,vce(cluster gwf_caseid)
		
		* Drop controls
		 xtset gwf_caseid year
		 qui xtprobit xonset lxyrs coldwar lt lnregion xpers2,vce(cluster gwf_caseid)	
		 est store c1
		 qui xtprobit xonset lxyrs xpers2,vce(cluster gwf_caseid)
		 lincom xpers2 
		 est store c2
		 qui xtprobit xonset lxyrs lt xpers2,vce(cluster gwf_caseid)
		 lincom xpers2 
		 est store c3
		 qui xtprobit xonset lxyrs lnregion xpers2,vce(cluster gwf_caseid)
		 lincom xpers2		 
		 est store c4
		 qui xtprobit xonset lxyrs coldwar xpers2,vce(cluster gwf_caseid)
		 lincom xpers2
		 est store c5
		 qui xtprobit xonset lxyrs coldwar lnregion xpers2,vce(cluster gwf_caseid)
		 lincom xpers2
		 est store c6
		 qui xtprobit xonset lxyrs coldwar lt lnregion xpers2,vce(cluster gwf_caseid)
		 lincom xpers2
		 est store c7
		 qui xtprobit xonset lxyrs lt lnregion xpers2,vce(cluster gwf_caseid)
		 lincom xpers2
		 est store c8	
		 qui xtprobit xonset lxyrs coldwar xpers2,vce(cluster gwf_caseid)
		 lincom xpers2
		 est store c9
		* No controls estimates table *
		estout c1 c2 c3 c4 c5 c6 c7 c8 c9 using TableB2.tex, cells(b(star  fmt(%9.3f)) se(par fmt(%9.2f))) ///
			stats(r2 N N_clust) style(tex) replace label starlevels(* 0.05) title(\label{tabB2}) 
		 
		**** Adding covariates, one at a time ****
		use temp-fe,clear
		spearman xpers2 v2x_clpol v2clkill v2cltort v2juhcind v2x_jucon v2x_frassoc_thick v2x_freexp_altinf v2x_partipdem v2x_libdem v2x_polyarchy
		gen n=_n
		gen beta=.
		gen hi=.
		gen lo=.
		gen hi90=.
		gen lo90=.
		gen varname=""
		xtset gwf_caseid year
		global cvar ="coldwar lxyrs lt lnregion"
		global varA= "lag_xongoing loggdp logoil support gwf_mil leadermil seizure_coup coup12 election ythbul4 wdipopurbmi wditrade"
		global varB= "polparcomp polpolcomp polexconst v2x_polyarchy v2x_libdem v2x_partipdem v2x_freexp_altinf v2x_frassoc_thick "
		global varC= "v2x_clpol v2x_jucon v2juhcind e_v2x_neopat priordem legcomp excluded monoethnic  multiethnic civwar grow"
		global varD= "v2cltort v2clkill lag_repress leadermil milethnic_homo nmc_logmilper nmc_logmilex effectivenumber debruin_cbcount debruin_ha_cbcount lpop"
		local var = "$varA $varB $varC $varD"
		local i =1
		foreach v of local var {
			di "`v'"
			qui xtset gwf_caseid
			qui xtprobit xonset $cvar `v' xpers2,vce(cluster gwf_caseid)	 
			qui nlcom _b[xpers2],post
			matrix beta =e(b)  
			local b = beta[1,1]
			qui replace beta=`b' if n==`i'
			matrix var = e(V) 
			local se =var[1,1]
			qui replace hi = `b' + sqrt(`se')*1.96 if n==`i'
			qui replace lo = `b' - sqrt(`se')*1.96 if n==`i'
			qui replace hi90 = `b' + sqrt(`se')*1.65 if n==`i'
			qui replace lo90 = `b' - sqrt(`se')*1.65 if n==`i'
			qui replace varname = "`v'" if n==`i'
			local i = `i' +1
		}
		label define varlab 1 "Ongoing protest campaign" 2 "GDP per capita (log)" 3 "Oil per capita (log)"  ///
			4 "Support party" 5 "Military junta" 6 "Military leader" 7 "Coup seizure" 8 "Recent coup" 9 "Election" ///
			10 "Youth bulge" 11 "Urban pop." 12 "Trade"   13 "Parcomp" 14 "Polcomp" 15 "Exec. constr." ///
			16 "V-Polyarchy" 17 "V-Liberal" 18 "V-Participation" 19 "V-Free express. info" ///
			20 "V-Free express. org" 21 "V-Civil lib." ///
			22 "V-Judical constr." 23 "V-Judicial indep." 24 "V-Neopatrimonialism" 25 "Prior democracy" ///
			26 "Leg. comp." 27 "Excluded pop" ///
			28 "Monoethnic party" 29 "Multiethnic party" 30 "Civil war" 31 "Growth"  32 "V-Torture" ///
			33 "V-Kill" 34 "Repression" 35 "Military leader" ///
			36 "Ethnic homo military" 37 "Military personnel (log)" 38 "Military exp. (log)" ///
			39 "Effective number"  40 "Counter-weights" 41 "Heavily armed cw" 42 "Population",replace
		label values n varlab
		twoway (scatter beta n if n<=42,mcol(blue)yline(-.1275468,lcol(red)lpat(dash_dot))) ///
			(rspike hi lo n if n<=42,lw(vthin)lcol(blue)) ///
			(rspike hi90 lo90 n if n<=42,lcol(blue)lw(medium)ytitle("{&beta}{sub:Security personalization}", ///
			size(large)height(-2)) ///
			xtitle(Added covariate)yline(0,lpat(dash)lcol(gs6))xlab(1(1)42,valuelabel angle(90))legend(off))
		graph export "$dir\baseline-covariates.pdf", as(pdf)   replace

		* Calendar time trends *
		xtset gwf_caseid year
		xtprobit xonset lxyrs lt lnregion xpers2 if year<1989,vce(cluster gwf_caseid)	
		est store time1 
		xtprobit xonset lxyrs lt lnregion xpers2 if year>=1989,vce(cluster gwf_caseid)	
		est store time2
		xtprobit xonset lxyrs time* lt lnregion xpers2,vce(cluster gwf_caseid)	
		est store time3
 		xtprobit xonset lxyrs period* lt lnregion xpers2,vce(cluster gwf_caseid)	
		est store time4
 		xtprobit xonset lxyrs d19* d20  lt lnregion xpers2,vce(cluster gwf_caseid)	
		est store time5		
		
		* Logits *
		xtset gwf_caseid year
		logit xonset $cvar xpers2,vce(cluster gwf_caseid)
		est store logit1
		xtlogit xonset $cvar xpers2,vce(cluster gwf_caseid) 
		est store logit2
		replace lpop = lpop*5
		xtlogit xonset $cvar xpers2 m_coldwar m_lxyrs m_lt m_xpers2  m_lnregion,vce(cluster gwf_caseid)
		est store logit3
		margins,dydx(xpers2)
		firthlogit xonset $cvar xpers2 m_coldwar mp_lxyrs mp_lt mp_xpers2 mp_lnregion mp_xonset
		est store logit4
		margins,dydx(xpers2)
		
		* Cook et al approach to unit intercepts placed in a probit *
		 probit xonset i.fe $cvar xpers2,vce(cluster gwf_caseid)
		 
		* Drop each region *
		use temp-fe,clear
		xtset gwf_caseid year
		local region = "cacar casia ceeurope easia meast nafrica samerica ssafrica weu"
		local i =1
		foreach v of local region {
			qui xtprobit xonset $cvar xpers2 if region~="`v'",re vce(cluster gwf_caseid)
			lincom xpers2
			est store region`i'
			local i =`i' +1
		 }
		 
		* Drop each decade *
		  qui xtprobit xonset $cvar xpers2 if year>=1960,re vce(cluster gwf_caseid)
		  lincom xpers2
		  est store decade1
		  local decade = "1960 1970 1980 1990 2000"
		  local i =2
		  foreach v of local decade {
			qui xtprobit xonset $cvar xpers2 if d`v'==0,re vce(cluster gwf_caseid)
			lincom xpers2
			est store decade`i'
			local i =`i' +1
		  }
	  
		  * Duration time *
		  xtprobit xonset $cvar xpers2,vce(cluster gwf_caseid)
		  xtprobit xonset coldwar xyrs* lt lnregion xpers2,vce(cluster gwf_caseid)
		  gen lnXxpers =lxyrs*xpers2
		  qui xtprobit xonset coldwar lxyrs lt lnregion xpers2 lnXxpers,vce(cluster gwf_caseid)
		  test lnXxpers
		  gen y1Xxpers = xpers2*xyrs
		  gen y2Xxpers = xpers2*xyrs2
		  gen y3Xxpers = xpers2*xyrs3
		  qui xtprobit xonset coldwar xyrs* y1X y2X y3X lt lnregion xpers2,vce(cluster gwf_caseid)
		  test y1 y2 y3
		  
		  * Alternative coding of DV *
		  use temp-fe,clear
		  gen xonsetNAV13 =  nav13_start==1 
		  gen xonsetNAV21 =  nav21_start==1 
		  gen xonsetMEC = nvcstart 
		  * This comes from the Ulfelder and Chenoweth 2015 data:  
		  * "Major Episodes of Contention (MEC) data set to identify the onset of thesenonviolent episodes"
		  gen xonsetREN = xonset*renavco 
		  recode xonsetREN (.=0)
		  tab xonset xonsetNAV13
		  tab xonset xonsetNAV21
		  tab xonset xonsetMEC
		  tab xonsetMEC xonsetNAV13
		  tab xonsetMEC xonsetNAV21
		  
		  xtset gwf_caseid year
		  btscs xonsetMEC year cow,gen(yrsMEC)
		  btscs xonsetNAV13 year cow,gen(yrsNAV13)
		  btscs xonsetNAV21 year cow,gen(yrsNAV21)
		  btscs xonsetREN year cow,gen(yrsREN)

		  sum yrs*
		  gen lyrsNAV13 =ln(yrsNAV13+1)
		  gen lyrsNAV21 =ln(yrsNAV21+1)
		  gen lyrsMEC =ln(yrsMEC+1)
		  gen lyrsREN =ln(yrsREN+1)

		  gen yrsMEC2 = yrsMEC^2
		  gen yrsMEC3 = yrsMEC^3
		  gen yrsNAV212 = yrsNAV21^2
		  gen yrsNAV213 = yrsNAV21^3
		  gen yrsNAV132 = yrsNAV13^2
		  gen yrsNAV133 = yrsNAV13^3
		  gen yrsREN2 = yrsREN^2
		  gen yrsREN3 = yrsREN^3

		  qui:xtprobit xonsetMEC lyrsMEC period* lnregion lt lpop xpers2,vce(cluster gwf_caseid)
		  lincom xpers2
		  est store alt1
		  qui:xtprobit xonsetMEC yrsMEC* period* lnregion lt xpers2,vce(cluster gwf_caseid)
		  lincom xpers2
		  qui:xtprobit xonsetNAV13 lyrsNAV13 period* lnregion lt xpers2,vce(cluster gwf_caseid)
		  lincom xpers2
		  est store alt2
		  qui:xtprobit xonsetNAV13 yrsNAV13* period* lnregion lt xpers2,vce(cluster gwf_caseid)
		  lincom xpers2
		  qui:xtprobit xonsetNAV21 lyrsNAV21 period* lnregion lt xpers2,vce(cluster gwf_caseid)
		  lincom xpers2
		  est store alt3
		  qui:xtprobit xonsetNAV21 yrsNAV21* period* lnregion lt xpers2,vce(cluster gwf_caseid)
		  lincom xpers2
		  qui:xtprobit xonsetREN lyrsREN period* lnregion lt xpers2,vce(cluster gwf_caseid)
		  lincom xpers2
		  est store alt4
		  qui:xtprobit xonsetREN yrsREN* period* lnregion lt xpers2,vce(cluster gwf_caseid)
		  lincom xpers2
		  qui:xtprobit xonset  lxyrs period* lnregion lt xpers2,vce(cluster gwf_caseid)
		  est store alt0
		  
		  
		  **********************************
		  ************ Plots ***************
		  **********************************  
		  	* Reported RE & CRE models coefficient plot *
			coefplot (b0, msym(S)) (b1, msym(d)) (b2, msym(t)) (b3, msym(O)) (b4, msym(oh)), ///
				title("Personalization and Non-violent protest", size(small)) ///
				drop(_cons xyr* lxyr* m_* mXy_* y_* ) ///
				order(gwf_personal xpers xpers1 xpers2) xline(0) grid(glcolor(gs13)) ///
				mfcolor(white) xlabel(-.8(.2).4,labsize(vsmall))levels(95 90)  ///
				legend(lab(3 "Probit")lab(6 "RE probit")lab(9 "CRE" "probit") ///
				lab(12 "Two-way CRE" "probit")lab(15 "Interactive CRE"  "probit")  ///
				size(vsmall)pos(7)ring(0)) ysize(3) xsize(3.5) saving(r1, replace) ///
				xtitle("Coefficient estimate", size(small) height(6)) ///
				note("90 (thick) and 95 (thin) percent confidence intervals", size(vsmall) pos(6))
				graph export "$dir\base-models.pdf", as(pdf) replace
		  	
			* Personalism variables models coefficient plot *
			coefplot (p1, msym(d)) (p2, msym(t)) (p3, msym(oh)) (p4, msym(plus)) (p5, msym(P)) (p6, msym(T)), ///
				title("Personalization and Non-violent protest", size(small)) drop(_cons xyr* lxyr* ) ///
				order(gwf_personal xpers xpers1 xpers2) xline(0) grid(glcolor(gs13)) mfcolor(white) xlabel(-.4(.1).3,labsize(vsmall))  levels(95 90)  ///
				legend(lab(3 "Regime type" "logit") lab(6 "Personalization" "logit") lab(9 "Personalization"  "RE logit") lab(12 "Party pers." "RE logit") ///
				lab(15  "Security pers." "RE logit") lab(18 "Both types" "personalism") size(vsmall)) ysize(3) xsize(3.5) saving(r1, replace)	xtitle("Coefficient estimate", size(small) height(6)) ///
				note("90 (thick) and 95 (thin) percent confidence intervals", size(vsmall) pos(6))
				graph export "$dir\personalism-models.pdf", as(pdf)   replace
			
				
			* Logit models coefficient plot *
			coefplot (logit1, msym(d)) (logit2, msym(t)) (logit4, msym(plus)), ///
				title("Personalization and Non-violent protest", size(small)) drop(_cons lxyr* m* coldwar ) ///
				order(xpers2 lt lpopl lnregion) xline(0) grid(glcolor(gs13)) mfcolor(white) xlabel(-.75(.25).75,labsize(vsmall))  levels(95 90)  ///
				legend(lab(3 "Logit") lab(6 "RE" "logit")   lab(9 "Firth" "logit") size(vsmall))ysize(3) xsize(3.5) ///
				saving(r1, replace)	xtitle("Coefficient estimate", size(small) height(6)) note("90 (thick) and 95 (thin) percent confidence intervals", size(vsmall) pos(6))
				graph export "$dir\logit-models.pdf", as(pdf)   replace

				
				* LPMs(WFE from Imai and Kim 2019 uses R script pasted below) *
					 use temp-fe,clear
					 keep if sample==1
					 * set the binary treatment variable at the median of security personalism *
					 qui sum xpers2 if sample==1,detail
					 local m=r(p50)
					 gen treat1=xpers2>=(`m'-.0001) 
					    * Note: using the IRT-2PL instead of mixed yields same treatment variable 
						qui sum irtpers2 if sample==1,detail
						local m=r(p50)
						di `m'
						gen treat2=irtpers2>=(`m'-.0001) 
						tab treat2 treat1
						drop treat2 
					 xtset gwf_caseid year
					 gen l1treat1=l.treat1
					 gen l2treat1=l2.treat1
					 gen l3treat1=l3.treat1

					 gen l1xonset=l.xonset
					 gen l2xonset=l2.xonset
					 gen l3xonset=l3.xonset
					 
					 * Bivariate *
					 ttest xonset,by(treat1)
					 ttest xongoing,by(treat1)

					 * Comparison margins estimate from RE probit *
					 global cvarlpm="lnregion lt lxyrs"
					 xtset gwf_caseid year
					 qui xtprobit xonset period* $cvarlpm treat1,  vce(cluster gwf_caseid)
					 margins,dydx(treat1)
							
					 * 2-way FE *
					 qui reghdfe xonset $cvarlpm treat1,a(gwf_caseid year)vce(cluster gwf_leaderid)
					 lincom treat1
					 est store lpm1

					 * Unit-specific quadratic time trend *
					 set matsize 10000
					 xtset gwf_caseid year
					 xi:qui reg xonset i.gwf_caseid*time i.gwf_caseid*time2 ///
						lxyrs lt lnregion treat1,vce(cluster gwf_leaderid)
					 lincom treat1
					 est store lpm2
					 
					 * Interactive FE *
					 qui regife xonset $cvarlpm treat1,a(gwf_caseid year)factor(gwf_caseid year,1) vce(cluster gwf_leaderid)
					 lincom treat1
					 est store lpm3
					 
					 * Past treatment *
					 reghdfe xonset $cvarlpm treat1 l1treat1 l2treat1,a(gwf_caseid year)cluster(gwf_leaderid)
					 est store lpm4
					 lincom treat1
					 test l1treat1 l2treat1  

					 * Past outcome *
					 reghdfe xonset $cvarlpm treat1 l1xonset l2xonset,a(gwf_caseid year)cluster(gwf_leaderid)
					 est store lpm5
					 lincom treat1
					 test l1xonset l2xonset  
					 
					 * Past outcome, low level *
					 reghdfe xonset $cvarlpm treat1 l1v2cademmob l2v2cademmob,a(gwf_caseid year)cluster(gwf_leaderid)
					 est store lpm6
					 * Check with 4 lags, per Hamilton 2018 RES
					 reghdfe xonset $cvarlpm treat1 l1v2cademmob l2v2cademmob l3v2cademmob,a(gwf_caseid year)cluster(gwf_leaderid)
					 reghdfe xonset $cvarlpm treat1 l1v2cademmob l2v2cademmob l3v2cademmob l4v2cademmob,a(gwf_caseid year)cluster(gwf_leaderid)

					 * Check with continous *
					 reghdfe xonset $cvarlpm xpers2 l1v2cademmob l2v2cademmob l3v2cademmob,a(gwf_caseid year)cluster(gwf_leaderid)
					 reghdfe v2cademmob $cvarlpm treat1 l1v2cademmob l2v2cademmob l3v2cademmob,a(gwf_caseid year)cluster(gwf_leaderid)
					
					 * Past outcome and past treatment *
						* First test whether t-2 and t-3 lags of outcome/treatment are related to outcome conditional on 1-year lags *
					 reghdfe xonset $cvarlpm treat1 l2xonset l3xonset l2treat1 l3treat1,a(gwf_caseid year) cluster(gwf_leaderid)
 					 test l2xonset l3xonset l2treat1 l3treat1
					 reghdfe xonset $cvarlpm treat1 l1treat l1xonset l2xonset l3xonset l2treat1 l3treat1,a(gwf_caseid year) cluster(gwf_leaderid)
					 test l2xonset l3xonset l2treat1 l3treat1
					 ivreghdfe xonset $cvarlpm treat1 (l1treat1 l1xonset= l2xonset l3xonset l2treat1 l3treat1),a(gwf_caseid year) cluster(gwf_leaderid)
					 lincom treat1
					 est store lpm7
					 
					 * Check with continous treatment *
					 ivreghdfe xonset $cvarlpm xpers2 (l1.xpers2 l1xonset= l2xonset l3xonset l2.xpers2 l3.xpers2), ///
						a(gwf_caseid year) cluster(gwf_leaderid)
					 
					* LPM table *
					estout lpm* using TableB1.tex, cells(b(star  fmt(%9.3f)) se(par fmt(%9.2f))) ///
						stats(r2 N N_clust widstat) style(tex) replace label starlevels(* 0.05) title(\label{tabB1})		
					saveold temp-fe1,replace version(12)

						/*  R script for weight-FE (WFE) estimator: 
						ls()
						setwd("C:/Users/jgw12/Dropbox/Research/Pers-NAVCO/Pers-RENAVCO/data")
						getwd()

						rm(list=ls())
						set.seed(23789)

						## install missing packages, and update if newer version available

						packageList<-c("car","plm","xtable","wfe","sandwich","stargazer")

						for(i in 1:length(packageList)){
							if (!require(packageList[i],character.only = TRUE)){
								install.packages(packageList[i], repos="http://lib.stat.cmu.edu/R/CRAN/")
							}
						}
						# update.packages(ask = FALSE, dependencies = c('Suggests'), oldPkgs=packageList, repos="http://lib.stat.cmu.edu/R/CRAN/")
						 
						library(car)
						library(foreign)
						library(xtable)
						library(wfe)
						library(sandwich)
						library(stats4)
						 

						rm(list=ls())
						set.seed(89270258)
						data.temp<-read.dta("temp-fe1.dta")

						# Unweighted unit FE + period effects #
						mod.W1<-wfe(xonset~treat1+lxyrs+lnregion+lt+
									  period1+period2+period3+period4+period5+
									  period6+period7+period8+period9+period10+
									  period11+period12,
									 data=data.temp,treat="treat1",
									 unit.index="gwf_caseid",time.index="year",
									 method="unit", qoi="ate",hetero.se=TRUE,
									 auto.se=TRUE,unweighted=TRUE)
						 summary(mod.W1)
						 
						 # Weighted unit FE + period effects #
						 mod.W2<-wfe(xonset~treat1+lxyrs+lnregion+lt+
									   period1+period2+period3+period4+period5+
									   period6+period7+period8+period9+period10+
									   period11+period12,
									  data=data.temp,treat="treat1",
									  unit.index="gwf_caseid",time.index="year", 
									 method="unit",qoi="ate",hetero.se=TRUE,
									  auto.se=TRUE,unweighted=FALSE)
						summary(mod.W2)
						 

						# Two-way FE UNweighted #
						mod.W3<-wfe(xonset~treat1+lxyrs+lnregion+lt,
									data=data.temp,treat="treat1",
									unit.index="gwf_caseid",
									time.index="year", method="unit",
									qoi="ate",hetero.se=TRUE,
									auto.se=TRUE,unweighted = TRUE,
									estimator = "did")
						summary(mod.W3)

						# Two-way FE Weighted #
						mod.W4<-wfe(xonset~treat1+lxyrs+lnregion+lt,
									data=data.temp,treat="treat1",
									unit.index="gwf_caseid",
									time.index="year", method="unit",
									qoi="ate",hetero.se=TRUE,
									auto.se=TRUE,unweighted =FALSE,
									estimator = "did")
						summary(mod.W4)

						*/
						
			* Semiparametric estimator *
			use temp-fe,clear
			xtset gwf_caseid year
			set seed 897243574
			xi:qui xtsemipar xonset i.year xpers2 lt lxyrs lnregion,nonpar(xpers2)cluster(gwf_leaderid)deg(3)gen(a b)spline
			qui gen invb= invlogit(b) /* rescale on 0,1 */
			qui sum invb if sample==1
			di r(mean)
			qui sum xonset   if sample==1 & invb~=.
			replace invb = invb-(.49943119-r(mean)) if sample==1 
			sum invb xonset if sample==1 & invb~=.
			qui gen inva=invlogit(a)
			qui sum inva if sample==1
			di r(mean)
			qui sum xonset if sample==1 & inva~=.
			replace inva = inva-(.4983158-r(mean)) if sample==1 
			twoway (hist xpers2,bin(50)bcol(gs12)yscale(range(0 3000)axis(2))yaxis(2) ///	
				ylab(none,axis(2))freq ytitle("",axis(2)))(lpolyci invb xpers2,bw(.5) xtit(Security personalism)ylab(.03(.01).06) yscale(alt) ///
				legend(off) ytitle(Probability of campaign onset)col(blue)xlab(-1.5(.5)1.5)yline(.0390294)) 
 			graph export "$dir\onset-semipar.pdf", as(pdf) replace
			drop a b invb inva
			
			* Time trend models coefficient plot *
			coefplot (time1, msym(d)) (time2, msym(t)) (time3, msym(oh)) (time4, msym(plus)) (time5, msym(P)), ///
				title("Adjusting for time trends", size(small)) drop(_cons xyr* lxyr* time* d* coldwar period*) ///
				order(xpers2) xline(0) grid(glcolor(gs13)) mfcolor(white) xlabel(-.5(.25).5,labsize(vsmall))  levels(95 90)  ///
				legend(lab(3 "1946-1988" "period") lab(6 "1989-2010" "period") lab(9 "Five-year"  "periods") lab(12 "Non-linear"  "time trend")  ///
				lab(15 "Decade"  "dummies") size(vsmall)) ysize(3) xsize(3.5) saving(r1, replace)	///
				xtitle("Coefficient estimate", size(small) height(6))note("90 (thick) and 95 (thin) percent confidence intervals", size(vsmall) pos(6))
				graph export "$dir\time-models.pdf", as(pdf) replace
						
			* Drop each region *
			  local region = "cacar casia ceeurope easia meast nafrica samerica ssafrica weu"
			  local i =1
			  foreach v of local region {
				qui xtprobit xonset lxyrs lt lnregion xpers2 if region~="`v'",re vce(cluster gwf_caseid)
				lincom xpers2
				est store region`i'
				local i =`i' +1
			  }
			  
			* Drop each decade *
			  qui xtprobit xonset lxyrs lt lnregion xpers2 if year>=1960,re vce(cluster gwf_caseid)
			  lincom xpers2
			  est store decade1
			  local decade = "1960 1970 1980 1990 2000"
			  local i =2
			  foreach v of local decade {
				qui xtprobit xonset lxyrs lt lnregion xpers2 if d`v'==0,re vce(cluster gwf_caseid)
				lincom xpers2
				est store decade`i'
				local i =`i' +1
			  }
			  
			label var xpers "Personalization"
			label var xpers1 `""Party        " "personalization""'
			label var xpers2 `""Security    " "personalization""'
			label var lpopl "Population (log)"
			label var lt `""Leader" "tenure (log)""'
			label var lnregion `""Region NVC" "onsets (log)""'
			label var gwf_pers `""Personalist" "regime    ""'
			  
			* Drop regions models coefficient plot *
			coefplot (region1, msym(d)) (region2, msym(t)) (region3, msym(oh)) (region4, msym(plus)) (region5, msym(P)) ///
				(region6, msym(d)) (region7, msym(t)) (region8, msym(oh)) (region9, msym(plus)) , ///
				title("Drop one region at a time", size(small)) drop(_cons pyr* time* d* coldwar) ///
				order(xpers2) xline(0) grid(glcolor(gs13)) mfcolor(white) xlabel(-.3(.15).3,labsize(vsmall))  levels(95 90)  ///
				legend(lab(3 "Central" "America") lab(6 "Central" "Asia") lab(9 "Central/East"  "Europe") lab(12 "East"  "Asia") ///
				lab(15 "Middle"  "East") lab(18 "North"  "Africa") lab(21 "South"  "America") lab(24 "Sub-Saharan"  "Africa") ///
				lab(27 "Western"  "Europe") size(vsmall)) ysize(3) xsize(3.5) saving(r1, replace)	///
				xtitle("Coefficient estimate", size(small) height(6))note("90 (thick) and 95 (thin) percent confidence intervals", ///
				size(vsmall) pos(6))
				graph export "$dir\drop-region-models.pdf", as(pdf) replace
				
			* Drop decade models coefficient plot *
			coefplot (decade1, msym(d)) (decade2, msym(t)) (decade3, msym(oh)) (decade4, msym(plus)) (decade5, msym(P)) (decade6, msym(d)), ///
				title("Drop each decade at a time", size(small)) drop(_cons pyr*) ///
				order(xpers2) xline(0) grid(glcolor(gs13)) mfcolor(white) xlabel(-.3(.15).3,labsize(vsmall))  levels(95 90)  ///
				legend(lab(3 "1940s" "1950s") lab(6 "1960s") lab(9 "1970s") lab(12 "1980s") ///
				lab(15 "1990s") lab(18 "2000s") size(vsmall)) ysize(3) xsize(3.5) saving(r1, replace) xtitle("Coefficient estimate", ///
				size(small) height(6))note("90 (thick) and 95 (thin) percent confidence intervals", size(vsmall) pos(6))
				graph export "$dir\drop-decade-models.pdf", as(pdf) replace
				
			* Alternative DV models coefficient plot *
			coefplot (alt0, msym(d)) (alt1, msym(t)) (alt2, msym(oh)) (alt3, msym(plus)) (alt4, msym(S)), ///
				title("Security personalization and Non-violent protest" "campaign onsets, alternative DV tests", ///
				size(small)) drop(_cons xyrs* yrs* lyrs* period*) order(xpers2 lt lpopl lnregion) ///
				xline(0) grid(glcolor(gs13)) mfcolor(white) xlabel(-.3(.15).3,labsize(vsmall)) levels(95 90) ///
				legend(lab(3 "Original") lab(6 "MEC") lab(9 "NAVCO 1.3") lab(12 "NAVCO 2.1") lab(15 "NEVER") size(vsmall)) ///
				ysize(3) xsize(3.5) xtitle("Coefficient estimate", size(small) height(6))  ///
				note("90 (thick) and 95 (thin) percent confidence intervals", size(vsmall) pos(6))
				graph export "$dir\alt-models.pdf", as(pdf)   replace
				
		*** Cox model with shared frailty ***
			use temp-fe,clear
			gen oxyrs = xyrs+1
			sum oxyrs
			stset oxyrs, fail(xonset)
			stcox coldwar lt lnregion xpers2 lpop,shared(gwf_caseid) nohr
			estat phtest,detail
			gen lnoxyrs = ln(oxyrs)
			gen stXlt = lnoxyrs*lt
			stcox lt coldwar lnregion xpers2 lpop stXlt,shared(gwf_caseid) nohr
			estat phtest,detail	
			centile xpers2,centile(10 90)
			centile lnregion,centile(10 90)
			stcurve, hazard at1(xpers2=-1.6) at2(xpers2=1.3) title(Loyalty) ylab(0(.1).4)  ///
				xtitle(Time since last protest onset,height(6)) saving(r1.gph,replace) ///
				legend(lab(1 "Low loyalty") lab(2 "High loyalty") pos(5) col(1) ring(0)) range(1 40)  
			stcurve, hazard  at1(lnregion=-.7) at2(lnregion=1.8) title(Regional protest) ylab(0(.1).4) ///
				xtitle(Time since last protest onset,height(6)) saving(r2.gph,replace) ///
				legend(lab(1 "Low regional protest") lab(2 "High regional protest") pos(5) col(1) ring(0)) range(1 40) 
			gr combine r1.gph r2.gph,  
			graph export "$dir\cox-models.pdf", as(pdf)   replace
			
		*** IVSLS test ***
			use temp-fe,clear
			** Leader effects **
			xtset gwf_caseid year
			qui xtprobit xonset coldwar lxyrs lt lnregion xpers2, re vce(cluster gwf_caseid)
			margins,dydx(xpers2)
			xtset gwf_leaderid year
			qui xtlogit xonset coldwar lxyrs lt lnregion xpers2, re vce(cluster gwf_leaderid)
			margins,dydx(xpers2)
			qui xtlogit xonset coldwar lxyrs lt lnregion xpers2, fe  
			margins,dydx(xpers2)
				* Leader-specific time trend *
			xi:qui reg xonset i.gwf_leaderid*time lxyrs lt lnregion xpers2,cluster(gwf_leaderid)
			lincom xpers
			
			** IV **
			egen rcount = count(year),by(gwf_caseid)
			keep if rcount>=2
			xtset gwf_caseid year
			gen milit2 = militrank^2
			xtreg xonset i.year lxyrs lt lnregion xpers2,re vce(cluster gwf_caseid)
			est store iv1
			xtivreg xonset i.year lxyrs lt lnregion (xpers2=militrank milit2), ///
				re vce(cluster gwf_caseid)  ec2sls 
			est store iv2
				qui: xtreg xpers2 i.year lxyrs lt lnregion   militrank milit2,vce(cluster gwf_caseid)
				test militrank=0
				test milit2=0,ac
			reghdfe xonset lxyrs lt lnregion xpers2,a(gwf_caseid year) cluster(gwf_leaderid) 
			est store iv3
			xi:ivreg2h xonset i.year lxyrs lt lnregion (xpers2=militrank milit2), ///
				fe cluster(gwf_leaderid) gmm2s partial(i.year)
			est store iv5
			
			 * Within IV-probit *
			ivprobit xonset $cvar (xpers2=militrank milit2) m_coldwar m_lxyrs  m_lt m_xpers2 m_lnregion,vce(cluster gwf_leaderid)
			margins,dydx(xpers2)predict(pr)
			
			 ** IV tests table **
 			estout iv1 iv2 iv3 iv5 using TableB3.tex, cells(b(star  fmt(%9.4f)) se(par fmt(%9.4f))) ///
				stats(N N_clust widstat) style(tex) replace label starlevels(* 0.05) title(\label{tabB3})
		
		*******************************
		*** RMSE model fit analysis ***
		*******************************
			use temp-fe,clear
			xtset cowcode year
			gen l1poldurable = l.poldurable
			gen l1v2x_libdem=l.v2x_libdem
			gen l1ythbul = l.ythbul
			gen l1v2cltort=l1.v2cltort
			gen l1v2clkill=l1.v2clkill
			recode debruin_affcount (5 4 3=3)
			global seed = 2892742
			global k=10 /* folds Chenoweth and Ufelder use 5 folds*/
			global y="xonset"
			global varA= "lag_xongoing l1v2cademmob loggdp logoil support lt leadermil coup12 election ythbul4 wdipopurbmi wditrade"
			global varB= "polparcomp polpolcomp polexconst l1v2x_polyarchy l1v2x_libdem l1v2x_partipdem l1v2x_freexp_altinf l1v2x_frassoc_thick "
			global varC= "l1v2x_clpol l1v2x_jucon l1v2juhcind e_v2x_neopat lnregion legcomp excluded monoethnic  multiethnic civwar grow"
			global varD= "l1v2cltort l1v2clkill lag_repress l1poldurable milethnic_homo lag_milpersonnel lag_milspend effectivenumber debruin_cbcount debruin_ha_cbcount lpop"
			local var = "$varA $varB $varC $varD"
			sum `var'
			qui reg $y xyrs*
			keep if e(sample)==1  
			sutex year $y  $varA $varB $varC $varD xpers2, minmax    	
				xtset gwf_caseid year
			capture program drop crossvalrmse
			program define crossvalrmse
				mat A=r(est)
				mat U = J(rowsof(A),1,1)
				mat c = U'*A
				mat cm = c/rowsof(A)
				mat list cm
			end
			
			gen base=.
			gen rmse=.
			gen n=_n
			gen varname = ""
			local var = " $varA $varB $varC $varD xpers2"	
			local i =1
			foreach v of local var {
				di "`v'"
				qui set seed $seed
				qui crossfold reghdfe $y xyrs* if `v'~=.,a(gwf_caseid coldwar)  k($k) 
				crossvalrmse
				local r = cm[1,1]
				qui replace base = `r'   if n==`i'
				qui set seed $seed
				qui crossfold reghdfe $y xyrs* `v',a(gwf_caseid coldwar)  k($k) 
				crossvalrmse
				local r = cm[1,1]
				qui replace rmse = `r'   if n==`i'
				qui replace varname = "`v'" if n==`i'
				local i = `i' +1
			}
		 
			* Change in RMSE from baseline *
			gen ch = (((rmse)-(base))/(rmse))*100
			* Plot change in RMSE *
			sort ch
			gen xn=_n
			replace n = xn
			drop xn
			list n varname ch if n<=43,clean noobs
		
				
				*** Condition of baseline specification covariates ***
				global varA= "lag_xongoing l1v2cademmob loggdp logoil support   leadermil coup12 election ythbul4 wdipopurbmi wditrade"
				global varB= "polparcomp polpolcomp polexconst l1v2x_polyarchy l1v2x_libdem l1v2x_partipdem l1v2x_freexp_altinf l1v2x_frassoc_thick "
				global varC= "l1v2x_clpol l1v2x_jucon l1v2juhcind e_v2x_neopat   legcomp excluded monoethnic  multiethnic civwar grow"
				global varD= "l1v2cltort l1v2clkill lag_repress l1poldurable milethnic_homo lag_milpersonnel lag_milspend effectivenumber debruin_cbcount debruin_ha_cbcount lpop"	
				local var = "$varA $varB $varC $varD xpers2"	
				local i =1
				foreach v of local var {
					qui set seed $seed
					qui crossfold reghdfe $y xyrs* lnregion lt if `v'~=.,a(gwf_caseid coldwar)  k($k) 
					qui crossvalrmse
					local r = cm[1,1]
					qui replace base = `r'   if n==`i'
					qui set seed $seed
					qui crossfold reghdfe $y xyrs* lnregion lt `v',a(gwf_caseid coldwar)  k($k) 
					qui crossvalrmse
					local r = cm[1,1]
					qui replace rmse = `r'   if n==`i'
					qui replace varname = "`v'" if n==`i'
					local i = `i' +1
				}
				replace ch = (((rmse)-(base))/(rmse))*100
				sort ch
				replace n =_n
				label define vlab 1 "Election" 2  "{bf:Security personalism}" 3 "Monoethnic party"  ///
					4 "Legislative competitiveness" 5 "Multiethnic party" 6 "V-Civil liberties" 7 "V-Free expression info" ///
					8 "Civil conflict" 9 "V-Free expression orgs" 10 "V-Polyarchy"	11 "Log GDPpc" 12 "Support party" ///
					13 "V-Participation"  14 "Ethnically homogenous military" 15 "V-Judical constraint"  ///
					16 "V-Kill" 17 "V-Judicial independence" 18 "Ongoing protest campaign" 19 "V-Liberal democracy" ///
					20 "V-Torture" 21 "Oil per capita (log)" 22 "Military leader" 23 "Excluded pop" 24 "V-Neopatrimonialism" ///
					25 "Recent coup" 26 "Economic growth" 27 "State repression" 28 "Military spending (log)" ///
					29 "Military personnel (log)" 30"Population" 31 "Parcomp" 32 "Polcomp"  33 "Executive constraint"   ///
					34 "Prior pro-democracy mobilization" 35 "Polity durable" 36 "Urban population" ///
					37 "Youth buldge"  38 "Trade"  39 "Counter-weights" 40 "Heavily armed counter-weights" ///
					41 "Effective # military organizations",replace
				label values n vlab,
				list varname ch if n<=41,clean noobs
				twoway (scatter n ch if n<=41,ysize(10)msym(P)mcol(blue)xtitle("Change in RMSE") ///
					ytitle("",height(0)) xlab(0 (2) 10)  ///
					ylab(1(1)41,valuelabel)legend(off)xline(0,lcol(red)))
				graph export "$dir\rsme-onset.pdf", as(pdf)   replace
				
			*********************************************************
			*** Cross-Validation with 10 runs of 5-folds 		  ***
			*** for 41 variables added separately to the baseline ***
			*********************************************************	
				use temp-fe,clear
				drop groups
					gen l1poldurable = l.poldurable
					gen l1polparcomp = l.polparcomp
					gen l1polpolcomp = l.polpolcomp
					gen l1polexconst = l.polexconst
					gen l1v2x_libdem=l.v2x_libdem
					gen l1ythbul4 = l.ythbul4
					gen l1v2cltort=l1.v2cltort
					gen l1v2clkill=l1.v2clkill
					gen l1wdipopurbmi = l1.wdipopurbmi
					gen l1wditrade = l1.wditrade
					gen l1loggdp = l1.loggdp
					gen l1lpop = l1.lpop
				global s = 97454
				global y="xonset"
				global x="lxyrs lt lnregion coldwar"
				global k="5"
				global varA= "xpers2 l1v2cademmob l1loggdp l1lpop election l1ythbul4 lag_milpersonnel lag_milspend l1wdipopurbmi l1wditrade "
				global varB= "lag_xongoing l1polparcomp l1polpolcomp l1polexconst l1v2x_polyarchy l1v2x_libdem"
				global varC= "l1v2x_partipdem l1v2x_freexp_altinf l1v2x_frassoc_thick l1v2x_clpol l1v2x_jucon l1v2juhcind"
				global varD= "e_v2x_neopat legcomp excluded monoethnic multiethnic civwar grow l1v2cltort l1v2clkill lag_repress"
				global varE= "l1poldurable milethnic_homo  effectivenumber debruin_cbcount debruin_ha_cbcount logoil support leadermil coup12"
				
				gen auc0=.
				gen auc1=.
				gen ratio=.
				gen varname=""
				gen n=_n
				local var = "$varA $varB $varC $varD $varE"	
				local i =1
				foreach v of local var {
				    di "`v'"
					qui cvauroc $y $x if `v'~=.,k($k)seed($s)
					local r = r(mean_auc)
					qui replace auc0 = `r' if n==`i'
					qui cvauroc $y $x `v' if `v'~=.,k($k)seed($s)
					local r = r(mean_auc)
					qui replace auc1 = `r' if n==`i'
					qui replace varname = "`v'" if n==`i'
					qui xtset gwf_caseid year
					qui xtsum `v'
					qui replace ratio = r(sd_w)/(r(sd_w)+r(sd_b)) if n==`i'
					local i = `i' +1
				}
				gen ch = (((auc1)-(auc0))/(auc0))*100
				sort ch
				replace n =_n
				browse n varname auc0 auc1 ch 
				twoway (scatter ch ratio) (lpoly ch ratio)
				
				replace varname ="{bf:Security Personalism}" if varname=="xpers2"
				replace varname ="State repression" if varname=="lag_repress"
				replace varname ="Election" if varname=="nelda_election"
				replace varname ="Polity durable" if varname=="l1poldurable"
				replace varname ="V-Polyarchy" if varname=="l1v2x_polyarchy"
				replace varname ="V-Civil liberties" if varname=="l1v2x_clpol"
				replace varname ="Intnl conflict" if varname=="prio_lconflict_int_inter"
				replace varname ="Civil conflict" if varname=="prio_lconflict_int_intra"
				replace varname ="Ethnic military" if varname=="milethnic_homo"
				replace varname ="GDPpc" if varname=="l1loggdp"
				replace varname ="Counter balance" if varname=="debruin_affcount"
				replace varname ="# Security orgs" if varname=="effectivenumber"
				replace varname ="Youth bulge" if varname=="l1ythbul4"
				replace varname ="Population" if varname=="l1lpop"
				replace varname ="Urban population" if varname=="l1wdipopurbmi"
				replace varname ="Support party" if varname=="support"
				replace varname ="Excluded pop" if varname=="excluded"
				replace varname ="Prior coup" if varname=="coup12"
				replace varname ="Oil per capita (log)" if varname=="logoil"
				replace varname ="Legislative competitiveness" if varname=="legcomp"
				replace varname = "Multiethnic party" if varname=="multiethnic"
				replace varname = "Monoethnic party" if varname=="monoethnic"
				replace varname = "Trade" if varname=="l1wditrade"
				replace varname = "Military personnel" if varname=="lag_milpersonnel"
				replace varname = "Military spending" if varname=="lag_milspend"
				replace varname = "Prior pro-democracy mobilization" if varname=="l1v2cademmob"
				replace varname = "Election" if varname=="election"
				replace varname = "V-Free expression info" if varname=="l1v2x_freexp_altinf"
				replace varname = "Military leader" if varname=="leadermil"
				replace varname = "V-Free expression orgs" if varname=="l1v2x_frassoc_thick"
				replace varname = "Counter-weights" if varname=="debruin_cbcount"
				replace varname = "Heavily armed counter-weights" if varname=="debruin_ha_cbcount"
				replace varname = "Economic growth" if varname=="grow"
				replace varname = "Polity parcomp" if varname=="l1polparcomp"
				replace varname = "Polity polcomp" if varname=="l1polpolcomp"
				replace varname = "Polity exconst" if varname=="l1polexconst"
				replace varname = "V-Neopatrimonialism" if varname=="e_v2x_neopat"
				replace varname = "V-Participation" if varname=="l1v2x_partipdem"
				replace varname = "Prior ongoing protest" if varname=="lag_xongoing"
				replace varname = "Civil war" if varname=="civwar"
				replace varname = "V-Judicial independence" if varname=="l1v2juhcind"
				replace varname = "V-Judical constraint" if varname=="l1v2x_jucon"
				replace varname = "State torture" if varname=="l1v2cltort"
				replace varname = "State killing" if varname=="l1v2clkill"
				replace varname = "V-Liberal democracy" if varname=="l1v2x_libdem"		 
				twoway (scatter n ch if n>40 & n<=41,mlab(varname)mlabpos(left)msym(P)mcol(blue)) ///
					(scatter n ch if n<=40,mlab(varname)xscale(range(-4 12.5)) ///
					msym(P)mcol(blue)xlab(-4(4)12)legend(off)ylab(1(1)41,  ///
					labsize(tiny)valuelabel)ytitle(Added variables,size(small)) xline(0,lcol(red))  ///
					xtitle(Percent change in baseline test AUC,size(small)) title("Protest campaign onset")   ///
					text(6 9 "{it:Baseline AUC: 0.614}", size(small)) ysize(8))
				graph export "$dir\crossval-onset.pdf", as(pdf)   replace
	
				* Show that military variables don't increase AUC substantially when population is in the specification *
				local var ="lag_milsp lag_milper"
				foreach v of local var {
					qui cvauroc $y $x lpop if `v'~=.,k($k)seed($s)
					local r1 = r(mean_auc)
					qui cvauroc $y $x `v' lpop if `v'~=.,k($k)seed($s)
					local r2 = r(mean_auc)
					local ch = (`r2'-`r1')/(`r1')
					local s ="  "
					di "`v'" 
					di `r1'  
					di `r2' 
					di `ch'
				}
			
			

				
	******************************************************************
	**************************  REPRESSION  **************************
	******************************************************************
		***** Differenced repression model *****
 		use temp-fe,clear
		xtset gwf_caseid year
		gen ongXpers=xongoing*xpers2
		gen onsXpers = xonset*xpers2
		xi:reghdfe d.repress d.xongoing d.xpers2,absorb(period* lt) cluster(gwf_caseid)
		est store rd1
		xi:reghdfe d.repress d.xongoing d.xpers2 d.ongXpers,absorb(period* lt) cluster(gwf_caseid)
		est store rd2
		centile xpers2 if e(sample)==1,centile(10 90)
		lincom d.xongoing + d.ongX*-1.6
		lincom d.xongoing + d.ongX*1.3
		xi:reghdfe d.repress d.xongoing d.coupA d.coupS d.election d.civwar d.xpers2 lxyrs,absorb(period* lt) cluster(gwf_caseid)
		est store rd3
		xi:reghdfe d.repress d.xongoing d.coupA d.coupS d.election d.civwar d.xpers2 d.ongXpers lxyrs,absorb(period* lt) cluster(gwf_caseid)
		est store rd4
		lincom d.xongoing + d.ongX*-1.6
		lincom d.xongoing + d.ongX*1.3

			** Check split sample **
			xi:qui reghdfe d.repress d.xongoing d.coupA d.coupS d.election d.civwar lxyrs if xpers2~=.,absorb(period* lt) cluster(gwf_caseid)
			lincom d.xongoing
			qui sum xpers2 if e(sample)==1,detail
			local m = r(p50)
			xi:qui reghdfe d.repress d.xongoing d.coupA d.coupS d.election d.civwar lxyrs if xpers2>=`m',absorb(period* lt) cluster(gwf_caseid)
			lincom d.xongoing	
			est store rd5
			xi:qui reghdfe d.repress d.xongoing d.coupA d.coupS d.election d.civwar lxyrs if xpers2<`m',absorb(period* lt) cluster(gwf_caseid)
			lincom d.xongoing
			est store rd6
			
		* Differenced estimates table *
		estout rd1 rd2 rd3 rd4 rd5 rd6 using TableC1.tex, cells(b(star  fmt(%9.4f)) se(par fmt(%9.3f))) ///
			stats(r2 N N_clust) style(tex) replace label starlevels(* 0.05) title(\label{tabC1})
			
		* Additional FEs *
		xi:reg d.repress d.xongoing d.coupA d.coupS d.election d.civwar d.xpers2 d.ongXpers lt, cluster(gwf_caseid)
		xi:reghdfe d.repress d.xongoing d.coupA d.coupS d.election d.civwar d.xpers2 d.ongXpers lt,absorb(year cow) cluster(gwf_caseid)
		xi:reghdfe d.repress d.xongoing d.coupA d.coupS d.election d.civwar d.xpers2 d.ongXpers lt,absorb(year gwf_caseid) cluster(gwf_caseid)
		xtset gwf_leaderid year
		xi:reghdfe d.repress d.xongoing d.coupA d.coupS d.election d.civwar d.xpers2 d.ongXpers lt,absorb(year) cluster(gwf_caseid)
		lincom d.xongoing + d.ongX*-1.6
		lincom d.xongoing + d.ongX*1.3

 		***** Repression in onset year & during protest campaigns *****
 		use temp-fe,clear
		xtset gwf_caseid
		egen mean_repress=mean(repress),by(gwf_caseid)
		egen mx_repress=mean(repress) if xongoing==0,by(gwf_caseid)
		egen m_repress=max(mx_),by(gwf_caseid)
		global x1 = "lpop lnregion lt lxyrs"
		centile xpers2 if xonset==1,centile(50)
		centile xpers2,centile(50)

		qui sum xpers2 if lpopl~=.,detail
		local m=r(p50)
		gen treat1=xpers2>=`m' 		 
		 ciplot lag1repress if xonset==1,by(treat) xtitle(Security personalization,size(small)height(3)) ///
			 xlab(2  "Low"  5  "High"  ) note("")  ///
			 ytitle("Lagged repression",height(2)) title("Repression{sub:t-1}") saving(r1.gph,replace)
		 ciplot lag2repress if xonset==1,by(treat) xtitle(Security personalization,size(small)height(3)) ///
			 xlab(2  "Low"  5  "High"  ) note("")  ///
			 ytitle("Lagged repression",height(2)) title("Repression{sub:t-2}") saving(r2.gph,replace)
		 ciplot lag3repress if xonset==1,by(treat) xtitle(Security personalization,size(small)height(3)) ///
			 xlab(2  "Low"  5  "High"  ) note("")  ///
			 ytitle("Lagged repression",height(2)) title("Repression{sub:t-3}") saving(r3.gph,replace)
		 ciplot lag_repress if xonset==1,by(treat) xtitle(Security personalization,size(small)height(3)) ///
			 xlab(2  "Low"  5  "High"  ) note("")  ///
			 ytitle("Lagged repression",height(2)) title("Repression{sub:t-1:t-3}") saving(r4.gph,replace)
		gr combine r1.gph r2.gph r3.gph r4.gph, ycommon col(2) iscale(.7) ysize(5) xsize(5)
		graph export "$dir\repression-ttests.pdf", as(pdf)   replace

		gen beta = .
		gen n =_n
		gen hi = .
		gen lo=.
		
		krls repress xpers2 $x1 period* lag1repress lag2repress lag3repress if xonset==1  
			mat o=e(Output)
			replace beta = o[1,1] if n==3
			replace hi = o[1,1]+1.95*o[1,2]  if n==3
			replace lo = o[1,1]-1.95*o[1,2]  if n==3
		krls repress xpers2 $x1 period* lag1repress lag2repress if xonset==1
			mat o=e(Output)
			replace beta = o[1,1] if n==2
			replace hi = o[1,1]+1.95*o[1,2]  if n==2
			replace lo = o[1,1]-1.95*o[1,2]  if n==2
		krls repress xpers2 $x1 period* lag1repress if xonset==1
			mat o=e(Output)
			replace beta = o[1,1] if n==1
			replace hi = o[1,1]+1.95*o[1,2]  if n==1
			replace lo = o[1,1]-1.95*o[1,2]  if n==1
		krls repress xpers2 $x1 period* mean_repress if xonset==1 /*condition on average regime repression level in all years */
			mat o=e(Output)
			replace beta = o[1,1] if n==4
			replace hi = o[1,1]+1.95*o[1,2]  if n==4
			replace lo = o[1,1]-1.95*o[1,2]  if n==4
		krls repress xpers2 $x1 period* m_repress if xonset==1    /*condition on average regime repression level in non-NVC years */
			mat o=e(Output)
			replace beta = o[1,1] if n==5
			replace hi = o[1,1]+1.95*o[1,2]  if n==5
			replace lo = o[1,1]-1.95*o[1,2]  if n==5
		krls repress xpers2 $x1 period* if xonset==1    			/*no condition on lagged repression */
			mat o=e(Output)
			replace beta = o[1,1] if n==6
			replace hi = o[1,1]+1.95*o[1,2]  if n==6
			replace lo = o[1,1]-1.95*o[1,2]  if n==6
		browse beta hi lo n if n<=6
		label define varlab 1 "t-1" 2 "t-1,t-2" 3 "t-1,t-2,t-3" 4 "Regime mean" 5 "All non-protest yrs" 6 "No lag",replace
		label values n varlab
		twoway (rspike hi lo n if n<=6,lcol(gs1)) (scatter beta n if n<=6,mcol(gs1)msym(O)yline(0,lcol(red)) yline(.0393769,lcol(blue)) ///
			ytitle("{&beta}{sub:Security personalization}",size(large))xtitle(Repression lag,height(3)) ///
			xscale(range(.75 6.25))yscale(range(0 .1))ylab(0(.05).2) xlab(1(1)6,valuelabel angle(45))legend(off))	
		graph export "$dir\repression-lags.pdf", as(pdf) replace
	
		**** Reported Models ****
		xtset cowcode year
		gen lag4repress=l.lag3repress
		krls repress period* lag1repress lag2repress lag3repress xpers2  if xonset==1, 
				est store rep1
		krls repress period* lag1repress lag2repress lag3repress xpers2 $x1 if xonset==1, 
				est store rep2
		krls repress period* lag1repress lag2repress lag3repress xpers2 lreduration $x1  if xongoing==1, d(dx) 
				est store rep3
		krls repress period* lag1repress lag2repress lag3repress lag4repress  xpers2 $x1 if xonset==1, 
			estout rep1 rep2 rep3 using Table1.tex, cells(b(star  fmt(%9.4f)) se(par fmt(%9.3f))) ///
		stats(r2 N N_clust) style(tex) replace label starlevels(* 0.05) title(\label{tab1})
		twoway (hist lreduration,col(gs14)yscale(range(0 20)axis(2))yaxis(2)bin(30)ylab(0(0)0,axis(2)) ytitle("",axis(2))) ///
			(lpolyci dx_xpers2 lreduration,bw(.75) lcol(blue*1.2)lpat(solid)col(blue*.25)legend(off) ///
			xtitle("Protest campaign duration (years, log scale)", ///
			size(small)) ytitle(Marginal effect of security personalism) yline(0,lcol(red))yscale(alt)yscale(alt axis(2)) ///
			xlab(0 "1" .69 "2"  1.099 "3"   1.61 "5" 2.079 "8" 2.48 "12")) 	 
		drop dx_* beta n hi lo
		graph export "$dir\pers-repression.pdf", as(pdf) replace
		
		* Check with OLS RE *
		xtset gwf_caseid year
		reg repress lag1repress lag2repress lag3repress period* xpers2 $x1 if xonset==1,cluster(gwf_caseid)  /* 6.3 percent marginal */
		xtreg repress lag1repress lag2repress lag3repress period* xpers2 $x1 if xonset==1,cluster(gwf_caseid)  /* 5.5 percent marginal */
		krls repress lag1repress lag2repress lag3repress period* xpers2 $x1 if xonset==1,   /* 4.8 percent marginal */
		
		* Additional covariates *
		gen beta = .
		gen n =_n
		gen var =""
		gen hi = .
		gen lo=.
		global vn = 25
		recode debruin_ha_cbcount (9 8 7 6 5 4 3 =3)
		local c="lnmembers territorial lag_xongoing civwar election coup coupA coupS coup12 excluded ethfrac milethnic_homo milethnic_het logoil loggdp xpers1 support nmc_logmilex nmc_logmilper  debruin_counterbalancing debruin_cbcount debruin_ha_cbcount effective v2juhcind v2x_jucon"
		local i = 1
		foreach v of local c {
			di "`v'"
			qui xi:krls repress xpers2 period* lag1repress lag2repress lag3repress  $x1 `v' if xonset==1
			mat o=e(Output)
			replace beta = o[1,1] if n==`i'
			replace var = "`v'" if n==`i'
			replace hi = o[1,1]+1.95*o[1,2]  if n==`i'
			replace lo = o[1,1]-1.95*o[1,2]  if n==`i'
			local i  = `i'+1
		}
		label define varlab 1 "Protest size" 2 `""Territorial" "campaign""' 3 `""Ongoing NV" "campaign""' 4 "Civil war" 5 "Election" ///
			6 "Coup attempt" 7 "Failed coup" 8 "Successful coup" 9 "Lag coup counts" 10 "Excluded pop." 11 "Ethnic fraction."   ///
			12 `""Ethnically" "homo military""' 13 `""Ethnically" "hetero military""' 14 "Oil per capita (log)"  15 "GDP per capita (log)"  ///
			16 "Party personalism" 17 "Support party" 18 "Mil. spending" 19 "Mil. personnel" 20  "Counterbalancing"  ///
			21  `""Counterbalance" "count""' 22  `""Counterbalance" "Heavily armed""' 23 `""Effective #" "mil. orgs""' ///
			24  "Judicial indep." 25  "Judicial constraint",replace
		label values n varlab
		twoway (rspike hi lo n if n<=$vn) (scatter beta n if n<=$vn,msym(O)mcol(gs1)yline(0,lcol(red)) yline(.0287,lcol(blue)) ///
			ytitle("{&beta}{sub:Security personalization}",size(large))xtitle(Added variable,height(-6)) ///
			xscale(range(.75 $vn.25))yscale(range(0 .08))ylab(0(.02).08) xlab(1(1)$vn,valuelabel angle(90))legend(off)xsize(8)ysize(4))		
		graph export "$dir\added-variables-repression.pdf", as(pdf) replace
		qui: reg repress lag_repress period* xpers2 $x1 if xonset==1,cluster(gwf_caseid)
		lincom xpers2
		qui: xtreg repress lag_repress period* xpers2 $x1 if xonset==1,cluster(gwf_caseid)
		lincom xpers2

		  * Check with VDem data on repression *
		  spearman vkill repress 
		  spearman vkill repress if xonset==1
		krls vkill period* lag_vkill xpers2 if xonset==1 
				est store rep1a
		krls vkill period* lag_vkill xpers2 $x1 if xonset==1
				est store rep2a
		krls vkill period* lag_vkill xpers2 lreduration $x1 if xongoing==1,d(dx)
				est store rep3a
		twoway lpolyci dx_xpers2 lreduration,bw(.7) legend(off) xtitle("Protest campaign duration (log)", ///
			size(small)) ytitle(Marginal effect of security personalism) yline(0,lcol(red))
		drop dx_*
		estout rep1a rep2a rep3a using TableC2.tex, cells(b(star  fmt(%9.4f)) se(par fmt(%9.3f))) ///
		stats(r2 N N_clust) style(tex) replace label starlevels(* 0.05) title(\label{tabC2})
		
		* Check with additive *
		xtset cow year
		gen ladd= l.additive
		krls additive period* ladd xpers2 $x1 if xonset==1

		**** Test Repression model with other forms of personalization *****  adjust Table 1, column 2 model 
		krls repress period* lag1repress lag2repress lag3repress gwf_pers $x1 if xonset==1
		est store p1
		krls repress period* lag1repress lag2repress lag3repress xpers $x1 if xonset==1
		est store p2
		krls repress period* lag1repress lag2repress lag3repress xpers1 $x1 if xonset==1
		est store p3
		krls repress period* lag1repress lag2repress lag3repress xpers2 $x1 if xonset==1
		est store p4
		krls repress period* lag1repress lag2repress lag3repress xpers1 xpers2 $x1 if xonset==1
		est store p5
		estout p1 p2 p3 p4 p5  using TableC3.tex, cells(b(star  fmt(%9.4f)) se(par fmt(%9.3f))) ///
		stats(r2 N N_clust) style(tex) replace label starlevels(* 0.05) title(\label{tabC3})
					
		*** Repression ECMs ***
		use temp-fe,clear
		gen ongXpers=xongoing*xpers2
		gen onsXpers = xonset*xpers2
		gen nvtime =lreduration
		recode nvtime (.=0)
		local var  = "repress xongoing ongXpers xpers2 civwar election coup coupS coupA lnregion lpop"
		foreach v of local var {
			xtset gwf_caseid year
			qui gen d`v'=d.`v'
			qui gen l`v'=l.`v'
		}
		
		* ECMs *
		global x= "dcivwar lcivwar delection lelection dcoupS lcoupS dcoupA lcoupA"
		reghdfe drepress lrepress dxongoing lxongoing dxpers2 lxpers2 lt lxyrs, ///
			absorb(period* nvtime) vce(cluster gwf_caseid nvtime)
		centile xpers2 if e(sample)==1,centile(50)
		est store r1
		reghdfe drepress lrepress dxongoing lxongoing dxpers2 lxpers2 lt lxyrs $x, ///
			absorb(period* nvtime) vce(cluster gwf_caseid nvtime)
		est store r2
		reghdfe drepress lrepress dxongoing dongXpers lxongoing longXpers dxpers2 lxpers2 lt lxyrs, ///
			absorb(period* nvtime) vce(cluster gwf_caseid nvtime)
		est store r3
		reghdfe drepress lrepress dxongoing dongXpers lxongoing longXpers dxpers2 lxpers2 lt lxyrs $x, ///
			absorb(period* nvtime) vce(cluster gwf_caseid nvtime)
		est store r4
		
				* Check ECM split sample *
				qui reghdfe drepress lrepress dxongoing lxongoing  lt lxyrs if xpers2>.0260124,  absorb(period* nvtime) vce(cluster gwf_caseid)
				lincom dxongoing
				qui reghdfe drepress lrepress dxongoing lxongoing  lt lxyrs if xpers2< .0260124,  absorb(period* nvtime) vce(cluster gwf_caseid)
				lincom dxongoing
				qui reghdfe drepress lrepress dxongoing lxongoing  lt lxyrs $x if xpers2> .0260124,  absorb(period* nvtime) vce(cluster gwf_caseid)
				lincom dxongoing
				qui reghdfe drepress lrepress dxongoing lxongoing  lt lxyrs $x if xpers2<.0260124,  absorb(period* nvtime) vce(cluster gwf_caseid)
				lincom dxongoing
				
		 * ECM plots *
			interflex drepress dxongoing xpers2 lrepress lxongoing lnregion lpop lt lxyrs, ///
				fe(period* nvtime) cluster(gwf_caseid) type(kernel)  bw(5)
			mat e=r(margeff)
			gen g=_n
			gen b=.
			gen se=.
			gen npers=.
 			forval i = 1(1)50 {
  				qui replace npers=e[`i',1] if g==`i'
  				qui replace b=e[`i',2] if g==`i'
  				qui replace se=e[`i',3] if g==`i'
			}
			gen hi5=.
			gen hi10=.
			gen lo5=.
			gen lo10=.
			replace hi5  = b+1.96*se
			replace hi10 = b+1.65*se
			replace lo5  = b-1.96*se
			replace lo10 = b-1.65*se
			graph twoway (hist xpers2,bin(50)freq bcolor(gs14)yscale(off)) ///
				(rspike hi5 lo5 npers if g<=50,lcolor(blue*.35)lwidth(thin)yaxis(2)yscale(alt axis(2)) ///
				xtitle("Security personalization",size(small)height(4)) ///
				yline(0,lp(dash)axis(2)) ytitle("Marginal effect of Protest campaign on Repression", ///
				size(small) height(4)axis(2))ylab(-.12(.04).12,axis(2)) ///
				title("Short-run",size(med))) ///
				(rspike hi10 lo10 npers if g<=50,lcolor(blue*.35)lwidth(medthick)yaxis(2)yscale(alt axis(2))) ///
				(line b npers if g<=50,yaxis(2)yscale(alt axis(2))lpattern(solid)lcolor(blue*1.1) ///
				xscale(range(0 1)) xlab(-1.6(.8)1.6)legend(lab(1 "Security personalization") ///
				lab(2 "95% CI") lab(4 "Marginal effect") size(small)pos(6)col(3)ring(1)) ///
				legend(order(1 2 4))saving(g1.gph,replace)) 
			drop se b g npers 
			interflex drepress lxongoing xpers2 lrepress dxongoing lnregion lpop lt lxyrs, ///
				fe(period* nvtime) cluster(gwf_caseid) type(kernel)  bw(5)
			mat e=r(margeff)
			gen g=_n
			gen b=.
			gen se=.
			gen npers=.
 			forval i = 1(1)50 {
  				qui replace npers=e[`i',1] if g==`i'
  				qui replace b=e[`i',2] if g==`i'
  				qui replace se=e[`i',3] if g==`i'
			}
			replace hi5  = b+1.96*se
			replace hi10 = b+1.65*se
			replace lo5  = b-1.96*se
			replace lo10 = b-1.65*se
			graph twoway (hist xpers2,bin(50)freq bcolor(gs14)yscale(off)) ///
				(rspike hi5 lo5 npers if g<=50,lcolor(blue*.35)lwidth(thin)yaxis(2)yscale(alt axis(2)) ///
				xtitle("Security personalization",size(small)height(4)) ///
				yline(0,lp(dash)axis(2)) ytitle("Marginal effect of Protest campaign on Repression", ///
				size(small) height(4)axis(2))ylab(-.12(.04).12,axis(2)) ///
				title("Long-run",size(med))) ///
				(rspike hi10 lo10 npers if g<=50,lcolor(blue*.35)lwidth(medthick)yaxis(2)yscale(alt axis(2))) ///
				(line b npers if g<=50,yaxis(2)yscale(alt axis(2))lpattern(solid)lcolor(blue*1.1) ///
				xscale(range(0 1)) xlab(-1.6(.8)1.6)legend(lab(1 "Security personalization") ///
				lab(2 "95% CI") lab(4 "Marginal effect") size(small)pos(6)col(3)ring(1)) ///
				legend(order(1 2 4))saving(g2.gph,replace)) 
		gr combine g1.gph g2.gph,xsize(8) 
		graph export "$dir\interflex-repression.pdf", as(pdf) replace  
		drop se b g npers 

		*** Onset year relative to post-onset years ***
		xi:krls repress xpers2 xonset civwar lpop lnmembers lt if xongoing==1,d(dx) 
		twoway lpolyci dx_xonset xpers2,legend(off) yline(0,lcol(red))
		drop dx*
			
		******* NAVCO repression **********
		use temp-fe,clear
		gen hipers = xpers2>0 if xpers2~=.
		tab any_repression if xonset==1

		 * Onset years *
		tab any_repression hipers if xonset==1,col
		krls any_repression xpers2 if xonset==1  
		krls any_repression xpers2 lnregion lpop lnmembers lt period* if xonset==1,d(k1)  
		krls any_repression xpers2 lnregion lpop lnmembers lt period* lag1repress lag2repress lag3repress if xonset==1
		twoway lpolyci k1_xpers2 year,bw(8) lcol(blue*1.2)lpat(solid)col(blue*.25)legend(off) ///
			xtitle("Year", size(small)) saving(h1.gph,replace) ///
			ytitle(Marginal effect of security personalism)ylab(0 .02 .04 .06) ///
			yline(0,lcol(red))  xlab(1950(10)2010)tit(Onset years only)
			
		* All campaign years *
		tab any_repression hipers if xongoing==1,col
		krls any_repression xpers2 if xongoing==1  
		krls any_repression xpers2 lreduration lnregion lpop lnmembers lt period* if xongoing==1,d(k2)  
		krls any_repression xpers2 lreduration lnregion lpop lnmembers lt period* lag1repress lag2repress lag3repress if xongoing==1
		twoway lpolyci k2_xpers2 lreduration,bw(.5) lcol(blue*1.2)lpat(solid)col(blue*.25)legend(off) ///
			xtitle("Protest campaign duration (years, log scale)", size(small)) saving(h2.gph,replace) ///
			ytitle(Marginal effect of security personalism)ylab(0 .02 .04 .06) tit(All campaign years) ///
			yline(0,lcol(red))  xlab(0 "1" .69 "2"  1.39 "4"    2.079 "8" 2.48 "12")
		gr combine h1.gph h2.gph,xsize(8)
		graph export "$dir\navco-repression.pdf", as(pdf)   replace
		
	******************************************************
	***** Security forces defection during campaign ******
	******************************************************
		use temp-fe,clear
		xtset gwf_caseid year
		qui sum lnmembers if any_sec_def~=.
		gen xmembers = (lnmembers-r(mean))/ r(sd)
		qui sum xpers2
		gen ipers2 = (xpers2+abs(r(min)))/(r(max)+abs(r(min)))
		global x1 = "lt period2-period12 lpop xmembers lnregion repress lreduration"
		xi:logit any_sec_defect  xpers2 period2-period12,  vce(cluster id)
		est store def1
		xi:logit any_sec_defect  xpers2 lt period2-period12,  vce(cluster id)
		est store def2
		xi:logit any_sec_defect  xpers2 $x1,  vce(cluster id)
		est store def3
		xi:logit any_sec_defect  xpers2 $x1 milethnic_homo,  vce(cluster id)
		margins,dydx(xpers2 xmembers milethnic_homo) /* lnmembers & xpers2 both now Stdev==1 scale for margins */
		est store def4
		xi:logit any_sec_defect  ipers2 $x1 milethnic_homo,  vce(cluster id)
		sum ipers2 xmembers milethnic_homo if e(sample)==1
		margins,dydx(ipers2 xmembers milethnic_homo)  /* milethnic_homo & xpers2 both now (0,1) scale for margins */
		krls any_sec_defect xpers2 
		xi:krls any_sec_defect xpers2 $x1
		local c="lag_xongoing election civwar priordem coup coupA coupS excluded milethnic_het nmc_logmilper nmc_logmilex logoil loggdp xpers1 support"
		foreach v of local c {
			di "`v'"
			qui xi:logit any_sec_defect xpers2 $x1 i.region `v', cluster(id)
			lincom xpers2
			qui xi:krls any_sec_defect xpers2 $x1 i.region `v'
			matrix b = e(b)
			di b[1,1]
		}
		  * modeling unit hetero increases variance but also increases estimate size; so NOT modeling unit effect yields LOWER estimates, which are reported *
		xi:logit any_sec_defect xpers2 $x1,cluster(gwf_caseid)
		xi:xtlogit any_sec_defect xpers2 $x1,vce(cluster gwf_caseid)
		xi:xtlogit any_sec_defect xpers2 lt coldwar lpop xmembers lnregion repress lreduration,fe	
		estout def1 def2 def3 def4 using TableE1.tex, cells(b(star  fmt(%9.3f)) se(par fmt(%9.3f))) ///
		stats(r2 N N_clust) style(tex) replace label starlevels(* 0.05) title(\label{tabE1})
		
		**** Alternative measures of personalization ****
		xi:logit any_sec_defect gwf_pers $x1,vce(cluster id)
		est store def_p1
		xi:logit any_sec_defect xpers $x1,vce(cluster id)
		est store def_p2
		xi:logit any_sec_defect xpers1 $x1,vce(cluster id)
		est store def_p3
		xi:logit any_sec_defect xpers2 $x1,vce(cluster id)
		est store def_p4
		xi:logit any_sec_defect xpers1 xpers2 $x1,vce(cluster id)
		est store def_p5


	********************************************************************** 
	******************     Democratization outcomes     ******************
	**********************************************************************
	use temp-fe,clear
	tab gdem gwf_case_fail
	forval i =1/12{
		egen m_period`i'=mean(period`i'),by(gwf_caseid)
	}
 	
	* show within transform equivalency *
	qui reghdfe gdem ld lnregion xpers2,a(gwf_caseid  year)cluster(gwf_caseid)
	lincom xpers2
	est store dem1
	xtsum gdem if e(sample)==1
	qui reghdfe gdem ld lnregion period* xpers2,a(gwf_caseid)cluster(gwf_caseid)
	lincom xpers2
	qui xtreg gdem ld lnregion period* xpers2,fe i(gwf_caseid)cluster(gwf_caseid)
	lincom xpers2
	qui reg gdem ld lnregion period* xpers2 m_ld m_lnregion m_coldwar m_xpers2 m_gdem m_period*,cluster(gwf_caseid)
	lincom xpers2
	qui reg gdem ld lnregion period* xpers2 m_ld m_lnregion m_coldwar m_xpers2 m_period*,cluster(gwf_caseid)
	lincom xpers2	
 
	  * CRE * Mundlak and Chamberlain within transformation, See Wooldridge 2002
	qui meprobit gdem ld lnregion period* xpers2 m_ld m_lnregion m_xpers2 m_period* || gwf_caseid: ,vce(cluster gwf_caseid)
	margins,dydx(xpers2) 
	est store dem2
	xtsum gdem if e(sample)==1
	qui probit gdem ld lnregion period* xpers2 m_ld m_lnregion m_xpers2  m_period*,cluster(gwf_caseid)
	margins,dydx(xpers2)
	
	* Add lagged trend in democracy *
	reghdfe gdem l1v2x_pol l2v2x_pol ld lnregion xpers2,a(gwf_caseid year)cluster(gwf_caseid)
	est store dem3
	xtsum gdem if e(sample)==1
		* Check with 3 & 4 year lags *
	reghdfe gdem l1v2x_pol l2v2x_pol l3v2x_pol l4v2x_pol ld lnregion xpers2,a(gwf_caseid year)cluster(gwf_caseid)
	egen m_l1v2x_pol=mean(l1v2x_pol) if sample==1,by(gwf_caseid)
	egen m_l2v2x_pol=mean(l2v2x_pol) if sample==1,by(gwf_caseid)
	meprobit gdem l1v2x_pol l2v2x_pol ld lnregion period* xpers2 ///
		m_ld m_lnregion m_xpers2 m_l1v2x_pol m_l2v2x_pol m_period* ///
		|| gwf_caseid: ,vce(cluster gwf_caseid)
	xtsum gdem if e(sample)==1
	est store dem4
	margins,dydx(xpers2) 
	
	* Interactive fixed effects *
	regife gdem ld lnregion xpers2 ,a(gwf_caseid year) factor(gwf_caseid year,1) vce(cluster gwf_caseid)
	est store dem5
	
	* Nonlinear case-specific time trends *
	xi:reg gdem i.gwf_caseid*time i.gwf_caseid*time2 ld lnregion xpers2 ,  vce(cluster gwf_leaderid)
	est store dem6
	
	estout dem1 dem2 dem3 dem4 dem5 using Table2.tex, cells(b(star  fmt(%9.4f)) se(par fmt(%9.3f))) ///
		stats(r2 N N_clust) style(tex) replace label starlevels(* 0.05) title(\label{tab2})
 
	
	* Interactive fixed effects, within: m_ are panel unit means; y_ are year means; and mXy_ are their interaction *
	meprobit gdem ld lnregion  xpers2 ///
		m_ld m_lnregion m_xpers2 y_ld y_lnregion y_xpers2  ///
		mXy_ld mXy_lnregion mXy_xpers2 ||  ///
		gwf_caseid: ,vce(cluster gwf_caseid)

	 * Check with IRT-2PL instead of generalized SFP *
	 reghdfe gdem ld lnregion irtpers,a(gwf_caseid year)cluster(gwf_caseid)
	 reghdfe gdem l1v2x_pol l2v2x_pol ld lnregion irtpers,a(gwf_caseid year)cluster(gwf_caseid)

	* Kernel estimator for Figure 4 *
	/* krls gdem ld lnregion coldwar xpers2 m_ld m_lnregion m_coldwar m_xpers2, ///
		d(k)lambda(2.4)ltolerance(4.5)sigma(8)
	replace k_xpers2=k_xpers2*100 */
	use temp-krls,clear
	twoway lpolyci k_xpers2 year if year>=1949, yline(0,lcol(red)lpat(dash))legend(off)ylab(-2(1)0) ///
		xlab(1950(10)2010)xtit(Year) ytit(Marginal effect of Security personalism) ///
		tit(Security personalism and democratic transition) ///
		note("Kernel regression, within transformation to proxy for regime-case fixed effects",pos(6)size(vsmall))
	graph export "$dir\krls-dem.pdf", as(pdf) replace  
	*save temp-krls,replace
	
		
	*** Other ways of modeling the time trend ***
	use temp-fe,clear
	forval i =1/12{
		egen m_period`i'=mean(period`i'),by(gwf_caseid)
	}
	qui reghdfe gdem ld lnregion coldwar xpers2,a(gwf_caseid)cluster(gwf_caseid)
	est store timedem1
	lincom xpers2
	qui reghdfe gdem ld lnregion d19* d20 xpers2,a(gwf_caseid)cluster(gwf_caseid)
	est store timedem2
	lincom xpers2
	qui reghdfe gdem ld lnregion period* xpers2,a(gwf_caseid)cluster(gwf_caseid)
	est store timedem3
	lincom xpers2
	qui reghdfe gdem ld lnregion time xpers2,a(gwf_caseid)cluster(gwf_caseid)
	est store timedem4
	lincom xpers2
	qui reghdfe gdem ld lnregion time time2 xpers2,a(gwf_caseid)cluster(gwf_caseid)
	est store timedem5
	lincom xpers2
	xtset gwf_caseid year
	regife gdem ld lnregion xpers2,a(gwf_caseid year)factor(gwf_caseid year,1)vce(cluster gwf_caseid)
	est store timedem6
		  * Nonlinear models *
	qui xtprobit gdem ld lnregion xpers2 m_ld m_lnregion m_xpers2 y_ld y_lnregion y_xpers2,vce(cluster gwf_caseid)
	est store timedem7
	 local var = "d1960 d1970 d1980 d1990 d2000 time time2"
	 foreach v of local var {
		egen m_`v' =mean(`v') if sample==1,by(gwf_caseid)
	 }
	qui xtprobit gdem ld lnregion xpers2 m_ld m_lnregion m_xpers2 d19* d20 m_d19* m_d20,vce(cluster gwf_caseid)
	est store timedem8
	qui xtprobit gdem ld lnregion xpers2 m_ld m_lnregion m_xpers2 period*,vce(cluster gwf_caseid)
	est store timedem9
	qui xtprobit gdem ld lnregion xpers2 m_ld m_lnregion m_xpers2 time m_time,vce(cluster gwf_caseid)
	est store timedem10
	qui xtprobit gdem ld lnregion xpers2 m_ld m_lnregion m_xpers2 time time2 m_time m_time2,vce(cluster gwf_caseid)
	est store timedem11
	qui xtprobit gdem ld lnregion xpers2 m_ld m_lnregion m_xpers2 ///
		y_ld y_lnregion y_xpers2 mXy_ld mXy_lnregion mXy_xpers2,vce(cluster gwf_caseid) 
	est store timedem12
	
	label var xpers2 " "
	coefplot (timedem1, msym(d))(timedem2, msym(t))(timedem3, msym(oh))(timedem4, msym(plus))(timedem5, msym(P)) ///
			(timedem6, msym(d)),title("Fixed effects linear probability model", size(med))  relocate(xpers2 = 1)  ///
			drop(lnregion coldwar _cons ld m_* y_* mXy_* d19* d20* period* time*) vertical xtit(Model specification,size(small)height(5)) ///
			order(xpers2) yline(0) grid(glcolor(gs13)) mfcolor(white) ylabel(-.04(.02)0,labsize(vsmall)) ///
			 levels(95 90)   xlab(.64 `""Year" "effects""' .78 `""Decade" "effects""' .925 `""Period" "effects""' 1.07 `""Linear" "time trend""' ///
			1.21 `""Non-linear" "time trend""' 1.35 `""Interactive" "fixed effects""',labsize(small) )  ///
			ysize(3) xsize(3.5) xscale(range(0.62 1.35))saving(g1, replace) ytit("Marginal effect estimate", size(small) ///
			height(4))legend(off)note("90 (thick) and 95 (thin) percent confidence intervals",size(small)ring(1) pos(6))
		
			 
		gen n=_n
		gen beta=.
		gen hi=.
		gen lo=.
		gen hi90=.
		gen lo90=.
		local i=1
		forval v =7/12 {
			di "`v'"
			est restore timedem`v'	 
			margins,dydx(xpers2)predict(pu0)  
			matrix beta =r(b)  
			local b = beta[1,1]
			qui replace beta=`b' if n==`i'
			margins,dydx(xpers2)predict(pu0)post
			matrix var = r(V) 
			local se =var[1,1]
			qui replace hi = `b' + sqrt(`se')*1.96 if n==`i'
			qui replace lo = `b' - sqrt(`se')*1.96 if n==`i'
			qui replace hi90 = `b' + sqrt(`se')*1.65 if n==`i'
			qui replace lo90 = `b' - sqrt(`se')*1.65 if n==`i'
			local i = `i' +1
		}
		graph twoway (rspike hi lo n if n<=6,lcolor(blue*.5)lwidth(thin) ///
			xtitle("Model specification",size(small)height(6)) ///
			yline(0,lp(dash)) ytitle("Marginal effect estimate", ///
			size(small) height(5))ylab(-.04(.02)0,)title("Correlated random effects probit",size(med))) ///
			(rspike hi90 lo90 n if n<=6,lcolor(blue*1) lwidth(medthick)) ///
			(scatter beta n if  n<=6,msym(dot)lpattern(solid)mcolor(blue*1.5) ///
			xscale(range(0.75 6.25))legend(lab(1 "95 ci")lab(3 "Estimate") size(small)pos(6)col(3)ring(1)) ///
			legend(order(3 1))saving(g2.gph,replace) xlab(1 `""Year" "effects""' 2 `""Decade" "effects""' ///
			3 `""Period" "effects""' 4 `""Linear" "time trend""' ///
			5 `""Non-linear" "time trend""' 6 `""Interactive" "fixed effects""'))
		gr combine g1.gph g2.gph,xsize(6)ysize(3)
					graph export "$dir\time-dem.pdf", as(pdf) replace  

			
	* Alternative Fixed effects *
	qui reghdfe gdem ld lnregion xpers2,a(gwf_caseid year)cluster(gwf_caseid)
	est store fedem1
	lincom xpers2
	qui reghdfe gdem ld lnregion xpers2,a(cowcode year)cluster(gwf_caseid)
	est store fedem2
	lincom xpers2
	qui reghdfe gdem lt lnregion xpers2,a(gwf_leaderid year)cluster(gwf_leaderid)
	est store fedem3
	lincom xpers2
	qui xi:xtreg gdem i.year ld lnregion xpers2,i(gwf_caseid)cluster(gwf_caseid)
	est store fedem4
	lincom xpers2
	qui xi:xtreg gdem i.year ld lnregion xpers2,i(cowcode)cluster(cowcode)
	est store fedem5
	lincom xpers2
	qui xi:xtreg gdem i.year lt lnregion xpers2,i(gwf_leaderid)cluster(gwf_leaderid)
	est store fedem6
	lincom xpers2
	
	coefplot (fedem1, msym(d))(fedem2, msym(t))(fedem3, msym(oh))(fedem4, msym(plus))(fedem5, msym(P)) ///
		(fedem6, msym(d)),title("Linear probability models", size(med))   ///
		drop(lnregion ld lt _Iyear_* _cons)   xtit(Model specification,size(small)height(5)) ///
		order(xpers2) xline(0) grid(glcolor(gs13)) mfcolor(white) xlabel(-.04(.02)0,labsize(vsmall)) ///
		levels(95 90)  ysize(3) xsize(3.5)   xtit("Marginal effect estimate", size(small) ///
		height(4))legend(lab(3 "Regime FE")lab(6 "Country FE")lab(9 "Leader FE") ////
		lab(12 "Regime RE")lab(15 "Country RE")lab(18 "Leader RE") ) ///
		note("90 (thick) and 95 (thin) percent confidence intervals; year effects in all specifications",size(vsmall)ring(1) pos(6))
	graph export "$dir\dem-effects.pdf", as(pdf) replace  
	
 
 
 	* All regime collapse *
	gen fail = gwf_case_fail
	recode fail (1=0) if gdem==1
 	qui reghdfe fail ld lnregion xpers2,a(gwf_caseid  year)cluster(gwf_caseid)
	lincom xpers2
	est store fail1
	xtsum fail if e(sample)==1
	xtset gwf_caseid year
	qui xtprobit fail ld lnregion period* xpers2 m_ld m_lnregion m_xpers2 m_period*, vce(cluster gwf_caseid)
	margins,dydx(xpers2) 
	lincom xpers2
	est store fail2
	xtsum fail if e(sample)==1
	estout fail1 fail2 using TableD1.tex, cells(b(star  fmt(%9.4f)) se(par fmt(%9.3f))) ///
		stats(r2 N N_clust) style(tex) replace label starlevels(* 0.05) title(\label{tabD1})
  
 	*** Cox duration model ***
	use temp-id,clear
	stset gwf_case_duration, failure(gdem) id(gwf_caseid)
	stcoxkm if gwf_case_duration<=100, by(treat) obs1opts(lcol(red)lpat(solid)mcol(red)msym(P)c(J)) ///
		obs2opts(lcol(blue)lpat(solid)mcol(blue)msym(P)c(J))  ///
		pred1opts(lcol(red)lpat(dash)mcol(red)msym(Oh)c(J)) ///
		pred2opts(lcol(blue)lpat(dash)mcol(blue)msym(Oh)c(J)) ///
		legend(pos(1)ring(0))ties(efron)xtit("Regime duration, years") ///
		note(Binary treatment is Security personalism greater than the sample median value,size(vsmall)pos(6))
		graph export "$dir\dem-coxkm.pdf", as(pdf)   replace
	stcox period* lnregion xpers2,vce(cluster gwf_caseid)nohr
	est store demcox1
		* Cox, within transform *
	forval i =1/12{
		egen m_period`i'=mean(period`i'),by(gwf_caseid)
	}
	stcox lnregion period* xpers2 m_lnregion m_period* m_xpers2,vce(cluster gwf_caseid)	nohr
	est store demcox2
		* Cox, with shared frailty by regime *
	gen ldXlnregion=ld*lnregion
	stcox period* lnregion xpers2 ldXlnregion,   shared(gwf_caseid)nohr
	est store demcox3
	estat phtest, rank detail
		* Cox, with strata by country *
	stcox period* lnregion xpers2, vce(cluster gwf_caseid)  strata(cowcode)	nohr
	estat phtest, rank detail
	est store demcox4
	stcox period* lnregion lt treat, vce(cluster gwf_caseid)  strata(cowcode)nohr
	estat phtest, rank detail
	stcox period* lnregion xpers2 loggdp lpop logoil, vce(cluster gwf_caseid)  strata(cowcode)nohr
	est store demcox5
	estat phtest, rank detail
	stcox period* lnregion xpers2 loggdp lpop logoil l1v2x_polyarchy l2v2x_polyarchy, vce(cluster gwf_caseid)   strata(cowcode)nohr	
	est store demcox6
	estout demcox* using TableD2.tex, cells(b(star  fmt(%9.4f)) se(par fmt(%9.3f))) ///
		stats(r2 N N_clust) style(tex) replace label starlevels(* 0.05) title(\label{tabD2})
  	
 
		**** Adding covariates, one at a time ****
		use temp-fe,clear
		recode debruin_affcount (5 4 3=3)
		spearman xpers2 v2x_clpol v2clkill v2cltort v2juhcind v2x_jucon v2x_frassoc_thick v2x_freexp_altinf v2x_partipdem v2x_libdem v2x_polyarchy
		gen n=_n
		gen beta=.
		gen hi=.
		gen lo=.
		gen hi90=.
		gen lo90=.
		gen varname=""
		xtset gwf_caseid year
		qui reghdfe gdem ld lnregion xpers2,a(gwf_caseid  year)cluster(gwf_caseid)
		nlcom _b[xpers2],post
		matrix beta =e(b)  
		global b = beta[1,1]	
		global num =40
		global cvar ="ld lnregion"
		global varA= "lpop loggdp logoil support gwf_mil leadermil seizure_coup coup12 election ythbul4 wdipopurbmi wditrade"
		global varB= "polparcomp polpolcomp polexconst l1v2x_polyarchy l2v2x_polyarchy l1v2x_partipdem l1v2x_freexp_altinf l1v2x_frassoc_thick "
		global varC= "l1v2x_clpol l1v2x_jucon l1v2juhcind l1e_v2x_neopat priordem legcomp inst excluded monoethnic multiethnic civwar grow"
		global varD= "milethnic_homo nmc_logmilper nmc_logmilex effectivenumber debruin_cbcount debruin_ha_cbcount debruin_affcount lag_xongoing"
		local var = "$varA $varB $varC $varD"
		local i =1
		foreach v of local var {
			di "`v'"
			qui xtset gwf_caseid
			qui reghdfe gdem `v' ld lnregion xpers2,a(gwf_caseid  year)cluster(gwf_caseid)
			qui nlcom _b[xpers2],post
			matrix beta =e(b)  
			local b = beta[1,1]
			qui replace beta=`b' if n==`i'
			matrix var = e(V) 
			local se =var[1,1]
			qui replace hi = `b' + sqrt(`se')*1.96 if n==`i'
			qui replace lo = `b' - sqrt(`se')*1.96 if n==`i'
			qui replace hi90 = `b' + sqrt(`se')*1.65 if n==`i'
			qui replace lo90 = `b' - sqrt(`se')*1.65 if n==`i'
			qui replace varname = "`v'" if n==`i'
			local i = `i' +1
		}
		label define varlab 1 "Population" 2 "GDP per capita (log)" 3 "Oil per capita (log)"  ///
			4 "Support party" 5 "Military junta" 6 "Military leader" 7 "Coup seizure" 8 "Recent coup" 9 "Election" ///
			10 "Youth bulge" 11 "Urban pop." 12 "Trade"   13 "Parcomp" 14 "Polcomp" 15 "Exec. constr." ///
			16 "V-Polyarchy,t-1" 17 "V-Polyarchy,t-2" 18 "V-Participation,t-1" 19 "V-Free express. info,t-1" ///
			20 "V-Free express. org,t-1" 21 "V-Civil lib.,t-1" 22 "V-Judical constr.,t-1" 23 "V-Judicial indep.,t-1" ///
			24 "V-Neopatrimonialism,t-1" 25 "Prior democracy" ///
			26 "Leg. comp." 27 "Institutions" 28 "Excluded pop" 29 "Monoethnic party" 30 "Multiethnic party" ///
			31 "Civil war" 32 "Growth" 33 "Ethnic homo military"  34 "Military personnel (log)" ///
			35 "Military exp. (log)" 36 "Effective number" 37 "Counter-weights" 38 "Heavily armed counter-weights" ///
			39 "Affiliated paramilitaries" 40 "Ongoing protests,t-1",replace
		label values n varlab
		twoway (scatter beta n if n<=$num,mcol(blue)yline($b,lcol(red)lpat(dash_dot))) ///
			(rspike hi lo n if n<=$num,lw(vthin)lcol(blue)) ///
			(rspike hi90 lo90 n if n<=$num,lcol(blue)lw(medium)ytitle("{&beta}{sub:Security personalization}", ///
			size(large)height(4)) ylab(-.06(.02)0) ///
			xtitle(Added covariate)yline(0,lpat(dash)lcol(gs6))xlab(1(1)$num,valuelabel angle(90))legend(off))
		graph export "$dir\dem-baseline-covariates.pdf", as(pdf)   replace
		
		*** IV-2SLS ***
		use temp-fe,clear
		xtset gwf_caseid year
		forval i =1/4 {
			gen l`i'vdem = l`i'.v2x_polyarchy
		}
		gen militrank2=militrank^2
		xi:ivreghdfe gdem ld lnregion (xpers2=militrank militrank2),absorb(gwf_caseid year)cluster(gwf_caseid) 
		est store demiv1
		xi:ivreghdfe gdem ld lnregion l1vdem l2vdem (xpers2=militrank militrank2),absorb(gwf_caseid year)cluster(gwf_caseid)   
		est store demiv2
		
		xi:ivreg2h gdem i.year  ld lnregion (xpers2=militrank militrank2),fe cluster(gwf_caseid) partial(i.year) gmm2s
		est store demiv3
		xi:ivreg2h gdem i.year l1vdem l2vdem  ld lnregion (xpers2=militrank militrank2),fe cluster(gwf_caseid) partial(i.year) gmm2s
		est store demiv4
		
		* Check with 3- and 4-lags *
		xi:ivreg2h gdem i.year l1vdem l2vdem l3vdem  ld lnregion (xpers2=militrank militrank2),fe cluster(gwf_caseid) partial(i.year) gmm2s
		xi:ivreg2h gdem i.year l1vdem l2vdem l3vdem l4vdem  ld lnregion (xpers2=militrank militrank2),fe cluster(gwf_caseid) partial(i.year) gmm2s
		
		estout demiv* using TableD3.tex, cells(b(star  fmt(%9.4f)) se(par fmt(%9.3f))) ///
			stats(r2 N N_clust widstat) style(tex) replace label starlevels(* 0.05) title(\label{tabD3})
			
			
		xi:reghdfe gdem lnregion l1vdem l2vdem l3vdem l4vdem xpers2,absorb(gwf_case_duration gwf_caseid year)cluster(gwf_caseid)
		xi:reghdfe gdem lnregion lpop loggdp logoil l1vdem l2vdem l3vdem l4vdem l.xpers2 xpers2,absorb(gwf_case_duration gwf_caseid year)cluster(gwf_caseid)
		xi:reghdfe gdem lnregion lpop loggdp logoil l1vdem l2vdem l3vdem l4vdem xpers2,absorb(gwf_case_duration gwf_caseid year)cluster(gwf_caseid)
	
	
	********* VDEM democracy score ************
	  reghdfe v2x_polyarchy ld lnregion xpers2,a(gwf_caseid year)cluster(gwf_leaderid)
	  egen minyr= min(year),by(gwf_caseid)
	  gen ovdem = l1v2x_poly if year==minyr
	  egen ivdem = max(ovdem),by(gwf_caseid)
	  reghdfe v2x_polyarchy ivdem ld lnregion xpers2,a(cowcode year)cluster(gwf_leaderid)

log close
		
 
		  
		****************** THE END *****************
		
